## Replication File for: Lily Tsai and Yiqng Xu. "Outspoken Insiders:
## Political Connections and Citizen Participation in Authoritarian
## China"

## Sensitivity Analysis (Figure 3)


rm(list=ls(all=TRUE))
source("vcovCluster.r") ## for clustered SE
library(foreign)
library(lmtest)
library(sandwich)

########################################
# Urban
########################################


d<-read.dta("match_urban.dta",convert.factors=F)
names(d)

Y<-"compgov"
tr<-"govconn"
Covars<-c("male","ageb2","ageb3","ageb4","govoff","hschool","college","ccp","hukou")
Covnames<-c("male","age31-40","age41-50","age>50","official","high school","college","ccp","Hukou")
#covarnames<-c("male","age","age-squared","govoff","edu")
locality<-"distrid"

data<-d
gammaset<-seq(0,.5,.01)
alphaset<-seq(0,.5,.01)
biasrecord<-matrix(NA,nrow=length(gammaset)*length(alphaset),ncol=3)
biasrecord_square<-matrix(NA,nrow=length(gammaset),ncol=length(alphaset))
position<-1
for (i in 1:length(gammaset)){
  for (j in 1:length(alphaset)){
    gamma=gammaset[i]  #define this as effect of Z on y
    alpha=alphaset[j] #define this as (x'x)^{-1}(x'z) for position corresponding to treatment
    #i.e. coef for treatment from regressing z on all x's.
    biasrecord[position,1]=gamma
    biasrecord[position,2]=alpha
    biasrecord[position,3]=gamma*alpha
    biasrecord_square[i,j]=gamma*alpha
    position=position+1
  }
}

## Convert from bias to effect esimate:

# Regression
(Xf<-as.formula(paste(Y," ~ ",paste(c(tr,Covars),collapse=" + ")," + as.factor(",locality,")",sep="")))
lmout=lm(Xf,data=data)
est<-lmout$coefficients[tr]
se<-coeftest(lmout,vcov = vcovCluster(lmout, cluster = data[,locality]))[tr,2] # clustered SE for X
effectrecord_square=est-biasrecord_square

# lowerbound: an effect would have to be larger than this to remain significantly non-zero.
(lowerboundsig<-0+1.96*se)

#Now get gammas and alphas for covars used for comparison:
covar_gammas<-lmout$coefficients[Covars]
covar_alphas<-replicate(length(Covars),NA)
for (i in 1:length(Covars)) {
  # regress Z on X
  Zf<-as.formula(paste(Covars[i]," ~ ",paste(c(tr,setdiff(Covars,Covars[i])),collapse=" + ")," + as.factor(",locality,")",sep=""))
  reg<-lm(Zf,data=data)
  covar_alphas[i]<-reg$coefficients[tr]
}

# Plot
fsize<-1.5
par(mar=c(5,5,2,2))
contour(xlim=c(-.01,.5),ylim=c(-.01,.5),x=gammaset,y=alphaset,z=effectrecord_square,nlevels=11,
        xlab="Effect of Confounder on the Probability of Filing Complaints",
        ylab="Effect of Being Politically Connected on Confounder",cex.lab=1*fsize)
lines(gammaset,lowerboundsig/gammaset, col="red",lwd=2,lty=2)
text(x=-0.09+median(gammaset),y=.05+lowerboundsig/median(gammaset),srt=-50
     ,labels="Reject null, p<0.05", col="red", cex=0.8*fsize)
points(abs(covar_gammas),abs(covar_alphas),pch=17,col=1)
text(x=abs(covar_gammas)+.00,y=abs(covar_alphas)+0.015,labels=Covnames,
     col=1,cex=.8*fsize)

########################################
# Rural
########################################


d<-read.dta("match_rural.dta",convert.factors=F)
names(d)

Y<-"compl"
tr<-"govconn"
Covars<-c("male","ageb2","ageb3","ageb4","leader","pschool","mschool","hschool","ccp")
Covnames<-c("male","age41-50","age51-60","age>60","cadre/official","primary school","middle school","high school","ccp")
locality<-"v_id"
data<-d

gammaset<-seq(0,.5,.01)
alphaset<-seq(0,.5,.01)
biasrecord<-matrix(NA,nrow=length(gammaset)*length(alphaset),ncol=3)
biasrecord_square<-matrix(NA,nrow=length(gammaset),ncol=length(alphaset))
position<-1
for (i in 1:length(gammaset)){
  for (j in 1:length(alphaset)){
    gamma=gammaset[i]  #define this as effect of Z on y
    alpha=alphaset[j] #define this as (x'x)^-1(x'z) for position corresponding to treatment
    #i.e. coef for treatment from regressing z on all x's.
    biasrecord[position,1]=gamma
    biasrecord[position,2]=alpha
    biasrecord[position,3]=gamma*alpha
    biasrecord_square[i,j]=gamma*alpha
    position=position+1
  }
}

## Convert from bias to effect esimate:

# Regression
(Xf<-as.formula(paste(Y," ~ ",paste(c(tr,Covars),collapse=" + ")," + as.factor(",locality,")",sep="")))
lmout=lm(Xf,data=data)
est<-lmout$coefficients[tr]
se<-coeftest(lmout,vcov = vcovCluster(lmout, cluster = data[,locality]))[tr,2] # clustered SE for X
effectrecord_square=est-biasrecord_square

# lowerbound: an effect would have to be larger than this to remain significantly non-zero.
(lowerboundsig<-0+1.96*se)

#Now get gammas and alphas for covars used for comparison:
covar_gammas<-lmout$coefficients[Covars]
covar_alphas<-replicate(length(Covars),NA)
for (i in 1:length(Covars)) {
  # regress Z on X
  Zf<-as.formula(paste(Covars[i]," ~ ",paste(c(tr,setdiff(Covars,Covars[i])),collapse=" + ")," + as.factor(",locality,")",sep=""))
  reg<-lm(Zf,data=data)
  covar_alphas[i]<-reg$coefficients[tr]
}

# Plot
fsize<-1.5
par(mar=c(5,5,2,2))
contour(xlim=c(-.01,.5),ylim=c(-.01,.5),x=gammaset,y=alphaset,z=effectrecord_square,nlevels=11,
        xlab="Effect of Confounder on the Probability of Filing Complaints",
        ylab="Effect of Being Politically Connected on Confounder",cex.lab=1*fsize)
lines(gammaset,lowerboundsig/gammaset, col="red",lwd=2,lty=2)
text(x=-0.01+median(gammaset),y=-.03+lowerboundsig/median(gammaset),srt=-50
     ,labels="Reject null, p<0.05", col="red", cex=.8*fsize)
points(abs(covar_gammas),abs(covar_alphas),pch=17,col=1)
text(x=abs(covar_gammas)+.00,y=abs(covar_alphas)+0.015,labels=Covnames,
     col=1,cex=.8*fsize)

